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ABSTRACT 

The development of a benchmark example for cyclic delamination growth 
prediction is presented and demonstrated for a commercial code. The example is based 
on a finite element model of a Double Cantilever Beam (DCB) specimen, which is 
independent of the analysis software used and allows the assessment of the 
delamination growth prediction capabilities in commercial finite element codes. First, 
the benchmark result was created for the specimen. Second, starting from an initially 
straight front, the delamination was allowed to grow under cyclic loading in a finite 
element model of a commercial code. The number of cycles to delamination onset and 
the number of cycles during stable delamination growth for each growth increment 
were obtained from the analysis. In general, good agreement between the results 
obtained from the growth analysis and the benchmark results could be achieved by 
selecting the appropriate input parameters. Overall, the results are encouraging but 
further assessment for mixed-mode delamination is required. 


INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common 
practice to characterize the onset and growth of delaminations. In order to predict 
delamination onset or growth, the calculated strain energy release rate components are 
compared to interlaminar fracture toughness properties measured over a range from 
pure mode I loading to pure mode II loading. 
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The virtual crack closure technique (VCCT) is widely used for computing energy 
release rates based on results from continuum (2D) and solid (3D) finite element (FE) 
analyses and to supply the mode separation required when using the mixed-mode 
fracture criterion [1,2]. The virtual crack closure technique was recently implemented 
into several commercial finite element codes. As new methods for analyzing 
composite delamination are incorporated into finite element codes, the need for 
comparison and benchmarking becomes important since each code requires specific 
input parameters unique to its implementation. 

An approach for assessing the delamination propagation capabilities in 
commercial finite element codes under static loading was recently presented and 
demonstrated for VCCT for ABAQUS R [3]. In this recent paper, benchmark results 
were created for full three-dimensional finite element models of the Double Cantilever 
Beam (DCB) and the Single Leg Bending (SLB) specimen. Then, starting from an 
initially straight front, the delamination was allowed to propagate in the finite element 
model. The load-displacement relationship and the total strain energy obtained from 
the propagation analysis results and the benchmark results were compared and good 
agreements could be achieved by selecting the appropriate input parameters. Overall, 
the results were encouraging but it was determined that further assessment for mixed- 
mode delamination is required [3], 

The objective of the present study was to create a benchmark example, 
independent of the analysis software used, which allows the assessment of the 
delamination fatigue growth prediction capabilities in commercial finite element 
codes. At the beginning, a benchmark example is created based on incremental finite 
element models of a DCB specimen. To avoid unnecessary complications, 
experimental anomalies such as fiber bridging were not addressed. First, a sample 
material and cyclic loading was selected. Second, the number of cycles to 
delamination onset, No, was calculated from the mode I fatigue delamination growth 
onset data of the material. Third, the number of cycles during stable delamination 
growth, ANg, was obtained incrementally from the material data for mode I fatigue 
delamination propagation by using growth increments of Aa= 0. 1 mm. Fourth, the total 
number of growth cycles, N G , was calculated by summing over the increments AN G . 
Fifth, the corresponding delamination length, a, was calculated by summing over the 
growth increments Aa. Finally, for the benchmark case where results for delamination 
onset and growth were combined, the delamination length, a, was calculated and 
plotted versus an increasing total number of load cycles Nt=Nd+Ng. After creating the 
benchmark, the approach was demonstrated for the commercial finite element code 
ABAQUS® . Starting from an initially straight front, the delamination was allowed to 
grow based on the algorithms implemented into the commercial finite element 
software. Input control parameters were varied to study the effect on the computed 
delamination increase during cyclic loading. It was assumed that delamination length 
increase during cyclic loading obtained from finite element analysis should closely 
match the growth shown in the benchmark example. The benchmark enabled the 
selection of the appropriate input parameters that yielded good agreement between 
the results obtained from the growth analysis and the benchmark results. Once the 
parameters have been identified, they may then be used with confidence to model 
delamination growth for more complex configurations. 



METHODOLOGY 


The methodology for delamination propagation, onset and growth [4] was applied 
to the DCB specimen (as shown in Figure 1) to create the benchmark example. Since 
the required material property input data were not readily available in the open 
literature for a single material, a fictitious set of properties was constructed for this 
benchmarking exercise. Individual properties for commonly used graphite/epoxy tape 
materials were obtained from the open literature to create this set to represent a typical 
graphite epoxy composite. The material properties are given in Table I. 

Static Fracture Toughness 

The mode I fracture toughness (mixed-mode ratio Gn/Gf=0) is generated 
experimentally using the DCB tests (as shown in Figure 1) [5], A fracture toughness 
Gi c = 0. 17 kJ/m 2 was used in this benchmarking exercise [6]. 

Fatigue Delamination Growth Onset 

The number of cycles to delamination onset, No, can be obtained from the 
delamination onset curve plotted in Figure 2 [7, 8]. The onset curve (solid green line) 
is a power law fit 


of the experimental data (open, green circles) obtained from a DCB test using the 
respective standard for delamination growth onset [7]. 

An example of an applied energy release G max = 0 . 1 362 kJ/ni 2 resulting in A/j=l 50 
cycles to delamination onset was included in Figure 2 (dashed red lines). The 
significance of the data will be discussed later. 
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dimensions 


B 25.0 mm 

2h 3.0 mm 

2L 150.0 mm 

a 0 30.5 mm 
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& max /2 0.67 mm 
6 min /2 0.067 mm 
R 0.1 


f 10.0 Hz 


£ 
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Figure 1. Double Cantilever Beam Specimen (DCB). 


TABLE I. MATERIAL PROPERTIES 


E\\ = 139.4 GPa 
v \2 = 0.30 
G\2 = 4.6 GPa 

Unidirectional Graphite/Epoxy Prepreg [3] 

£22= 10.16 GPa £33 = 10.16 GPa 

vi3=0.30 V 23 = 0.436 

G 13 = 4.6 GPa G 23 = 3.54 GPa 

Gj c = 0.17 kJ/m 2 

Fracture Toughness Data [3] 

Guc = G/i/c =0.49 kJ/m 2 rj= 1 .62 

m o=0. 2023 

Delamination Growth Onset Data [ 8 ] - Figure 2 
w/=-0. 078924 

Delamination Growth Rate Data (Paris Law) [ 6 ] - Figure 3 

Gi c = 0.17 kJ/m 2 

Gfh = 0.06 kJ/m“ 

c=2.44 10 6 

«=10.61 



Figure 2. Delamination growth onset for DCB specimen. 


Fatigue Delamination Growth 

The number of cycles during stable delamination growth, Ng, can be obtained 
from the fatigue delamination propagation relationship (Paris Law) plotted in Figure 3 
[6], The delamination growth rate (solid purple line) can be expressed as a power law 
function 


da 

~dN 


G il 


( 2 ) 


where da/dN is the increase in delamination length per cycle and G max is the maximum 
energy release rate at the front at peak loading. The factor c and exponent n are 
obtained by fitting the curve to the experimental data (open black circles) obtained 
from DCB tests [6]. The critical energy release rate or fracture toughness, G/ c , was 
included in the plot of Figure 3 (blue solid vertical line). Since composites do not 
exhibit the same threshold behavior commonly observed in metals, a cutoff value, G t h, 
was chosen below which delamination growth was assumed to stop (green solid 
vertical line) [6]. This benchmarking exercise ignores branching or fiber bridging and 
hence the Paris Law was not normalized with the static R-curve. 

An example of an applied energy release G max = 0.8 G/ c was included in Figure 3 
(red square). The significance of the data and the red arrow will be discussed later. 



Figure 3. Delamination growth rate (Paris Law). 


SPECIMEN AND FATIGUE TEST DESCRIPTION 


For the current numerical investigation, the Double Cantilever Beam (DCB) 
specimen, as shown in Figure 1, was chosen since it is simple, only exhibits the mode 
I opening fracture mode and had been used previously to develop an approach to 
assess the quasi-static delamination propagation simulation capabilities in commercial 
finite element codes [3], To avoid unnecessary complications, experimental anomalies 
such as fiber bridging were not addressed. For the current study, a DCB specimen 
made of graphite/epoxy with an unidirectional layup, [0] 2 4, was modeled. The material 
properties are given in Table I. The material, layup, overall specimen dimensions 
including initial crack length, a, were identical to the specimen used earlier [3], 

For the cyclic loading of the specimen, guidance was taken from a draft standard, 
where it is recommended to start the test at a maximum displacement, d max , which 
causes the energy release rate at the front, G/ max , to reach initially about 80% of Gi c . 
The details are discussed in reference 6. 

= 0.8 (3) 

G IC 

The maximum load, P max , and maximum displacement, d max /2, were calculated 
using the known quadratic relationship between energy release rate and applied load 
or displacement 



6„,-8 m ,V08 (5) 


where P cr u and d crit are the critical values. For the current study, a critical energy 
release rate G/ c =0.17 kJ/m 2 was used and the critical values P crit and d crit (grey dashed 
lines) were obtained from the benchmark for static delamination propagation [3] 
shown in the load-displacement plot in Figure 4. The calculated maximum load, P max , 
and calculated maximum displacement, d max /2, are shown in Figure 4 (dashed red line) 
in relationship to the static benchmark case (solid grey circles and dashed grey line) 
mentioned above. During constant amplitude cyclic loading of a DCB specimen under 
displacement control, the applied maximum displacement, d max /2=0.61 mm, is kept 
constant while the load drops as the delamination length increases (solid red circles 
and solid red line). The energy release rate corresponding to an applied maximum 
displacement b„ MA -/2=0.67 mm was calculated for different delamination lengths a 
using equation (5). The energy release rate decreases with increasing delamination 
length, a, as shown in Figure 5 (solid red circles and solid red line). Delamination 
growth was assumed to stop once the calculated energy release rate drops below the 
cutoff value, G>/„ (green solid horizontal line). The static benchmark case (solid grey 
circles and dashed grey line in Figure 5), where the delamination propagates at 
constant Gi c (solid blue line in Figure 5) was included for comparison. 



load P, 


G, kJ/m 



applied opening displacement <5/2, mm 


Figure 4. Critical load-displacement behavior for a DCB specimen. 



Figure 5. Energy release rate - delamination length behavior for DCB specimen. 


A load ratio Z?=0.1 was selected for this investigation as discussed in detail in 
reference 6. The corresponding minimum load, P m i„, and minimum displacement, 
d m inl2, were calculated 


p 6 

R _ — mm_ _ — fmn_ =Q | ^ 5 = 0 . 1-6 ( 6 ) 

j-j e min max V/ 

max ^ max 

Further, a frequency /=10 Hz was selected for this investigation as discussed in detail 
in reference 6. 

In the analysis, the displacement <5/2 applied to each arm of the specimen is 
represented as a function of time, t 

5 / 2 = [a 0 + b x • sin co(t - 1 0 )] • d max / 2 (7) 

where 6 max /2=0.61 mm is the maximum displacement. The constants ao= 0.55, 
b/=0.45, the circular frequency oo=20jr=62.832 and the starting time hrO.025 are 
calculated from load ratio i?=0.1 and the frequency /=10 Hz for testing. 


CREATING A BENCHMARK EXAMPLE FOR GROWTH PREDICTION 
Fatigue Delamination Growth Onset 

The number of cycles to delamination onset, No, may be obtained by solving 
equation (1) for No- 


G = 


m 0 ■ N 


m t 

D 



N D = c l -G c 
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where c, = — . Values for the constants c/ and C 2 are shown in Figure 2. 
m x 

At the beginning of the test, the specimen is loaded initially so that the energy 
release rate at the front, Gi max , reaches about 80% of Gu corresponding to 
Gi max =(). 1 362 kJ/m 2 as shown in Figure 2 (horizontal dashed red line). From the 
delamination onset curve, the number of cycles to delamination onset is determined as 
Nd= 150 cycles (vertical dashed red line). 

Fatigue Delamination Growth 

The number of cycles during stable delamination growth can be obtained by 
solving equation (2) for Ng 

N e = JdN.f^GZda 


( 9 ) 



As mentioned above, the specimen is loaded initially so that the energy release rate 
at the front, G/ max , reaches about 80% of Gu corresponding to Gi max = 0.1362 kJ/m 2 in 
the current study as shown in the Paris Law plot of Figure 3 (solid red square). 

For practical applications, equation (2) can be replaced by an incremental 
equivalent expression 


A a 
AN 


G n 

„ 


( 10 ) 


where for the current study, increments of Aa= 0.1 mm were chosen. Starting at the 
initial delamination length ao= 30.5 mm, the energy release rates G, max were obtained 
for each increment, i, from the curve fit (solid red circles and solid red line) plotted in 
Figure 5. These energy release rate values were then used to obtain the increase in 
delamination length per cycle or growth rate Aa/AN from the Paris Law in Figure 3. 
The growth rate rapidly decreases with increasing delamination length a as indicated 
by the red arrow in Figure 3. The number of cycles during stable delamination growth, 
Ng, was calculated by summing the increments AN, 

No-2 AN '-'l- G ^' ia < u > 

i = 1 /=1 C 

where k is the number of increments. The corresponding delamination length, a, was 
calculated by adding the incremental lengths Aa to the initial length ao. 

k 

a = a 0 Aa = a 0 + k ■ Aa (12) 

i=i 


Combined Fatigue Delamination Onset and Growth 

For the combined case of delamination onset and growth, the total life, Nt, may be 
expressed as 


N t =N d +N c (13) 

where, N D , is the number of cycles to delamination onset and Ng, is the number of 
cycles during delamination growth [9]. For this combined case, the delamination 
length, a, is plotted in Figure 6 for an increasing number of load cycles Nt. For the 
first Ad= 1 50 cycles, the delamination length remains constant (horizontal red line), 
followed by a growth section where - over Ng cycles - the delamination length 
increases following the Paris Law (crosses and solid red line). Once a delamination 
length is reached where the energy release rate drops below the assumed cutoff value, 
Go,, (as shown in the Paris Law in Figure 3) the delamination growth no longer 
follows the Paris Law (dashed grey line) and stops (horizontal solid red line). 



delamination 
length 
a, mm 



Figure 6. Delamination onset and growth behavior for DCB specimen. 


FINITE ELEMENT MODELING 

A typical two-dimensional finite element model of a Double Cantilever Beam 
(DCB) specimen is shown in Figure 7a. The specimen was modeled with solid plane 
strain elements (CPE4) and solid plane stress elements (CPS4) in ABAQUS® 
Standard 6.8, 6.9 and 6.9EF. Along the length, all models were divided into different 
sections with different mesh refinement. The DCB specimen was modeled with six 
elements through the specimen thickness (2h) as shown in the detail of Figure 7b. The 
resulting element length at the delamination tip was Aa=0.5 mm. A finer mesh, 
resulting in Aa=0.25 mm, was also generated. Additionally, three coarser meshes with 
a reduced number of elements in the length direction were also generated resulting in 
Aa=1.0 mm, Aa=1.25 mm and Aa=1.67 mm. The models are shown in reference 6. 

The plane of delamination was modeled as a discrete discontinuity in the center of 
the specimen. For the analysis with ABAQUS 8 6.8, 6.9 and 6.9EF, the models were 
created as separate meshes for the upper and lower part of the specimens with 
identical nodal point coordinates in the plane of delamination [10]. Two surfaces (top 
and bottom surface) were defined to identify the contact area in the plane of 
delamination as shown in Figure 7b. Additionally, a node set was created to define the 
intact (bonded nodes) region. 

A typical three-dimensional finite element model of the DCB specimen is shown 
in Figure 7c. Along the length, all models were divided into different sections with 
different mesh refinement. A refined mesh was used in the center of the DCB 
specimen. Across the width, a unifonn mesh was used to avoid potential problems at 
the transition between a coarse and finer mesh [3]. Through the specimen thickness 
(. 2h ), six elements were used. The resulting element length at the delamination tip was 


Aa=0.5 mm. The specimen was modeled with solid brick elements (C3D8I) which had 
yielded excellent results in a previous studies [3]. Two coarser meshes with a reduced 
number of elements in the width and length directions, resulting in Aa=1.0 mm and 
Aa=2.0 mm, were also generated as shown in reference 6. 

Three models of the DCB specimen were also generated with continuum shell 
elements (SC8R) resulting in an element length at the delamination tip Aa=0.5 mm, 
Aa=1.0mm and Aa=2.0mm, respectively. In the z-direction, only one element was 
used to model the thickness of the specimen. These less-refined models were used to 
study the effect on performance (CPU time), computed load/displacement behavior 
and growth prediction in comparison with the more refined models discussed above. 
Details are discussed in reference 6. 

For all the analyses performed, the low-cycle fatigue analysis in ABAQUS® 
Standard 6.8, 6.9 and 6.9EF was used to model delamination growth at the interfaces 
in laminated composites [10], A direct cyclic approach is part of the implementation 
and provides a computationally effective modeling technique to obtain the stabilized 
response of a structure subjected to constant amplitude cyclic loading. Delamination 
onset and growth predictions are based on the calculation of the strain energy release 
rate at the delamination front using VCCT. To determine propagation, computed 
energy release rates are compared to the input data for onset and growth from 
experiments as discussed in the methodology section. During the analysis, at least one 
element length at the crack tip is released along the interface after each stabilized cycle 
[10]. 

For all analyses, the elastic constants, the input to define the fracture criterion, and 
the parameters for delamination onset and delamination growth (Paris Law) were 
kept constant. The elastic constants, the fracture toughness values and the parameters 
required for delamination onset and growth are given in Table I. The parameters to 
define the load frequency (/= 10 Hz), the load ratio (R=0. 1 ) as well as the minimum 
and maximum applied displacement (<5,„„/2=0.067 mm and d max /2=0.67 mm) were 
also kept constant during all analyses. To study the effect on the computed onset and 
growth behavior during the analysis: 

• The number of terms used to define a Fourier series was varied. A Fourier series 
is used during the execution of the ABAQUS* Standard to approximate the 
periodic cyclic loading. 

• The size of the initial time increment used in the analysis was varied. 

• The input required to define the cyclic loading was altered. 

• The release tolerance was varied. Once a user specified release tolerance 
(( G-G c )/G c > release tolerance ) is exceeded during the analysis in ABAQUS® 
Standard, a cutback operation is performed which reduces the time increment. The 
cutback reduces the degree of overshoot and improves the accuracy of the local 
solution. 

• The solution controls were varied. 

It was assumed that the computed onset and growth behavior should closely match 
the benchmark results established earlier. Setting the value of the input parameters 
correctly is often an iterative procedure, which will be discussed later. Further 
details about the required input parameters are discussed in reference 6 where sample 
input fdes are also provided. 



initial delamination front 



bonded nodes 


a. Deformed two-dimensional FE-model of a DCB specimen with initial delamination 




c. Deformed three-dimensional FE-model of a DCB specimen with initial delamination front 

Figure 7. Typical finite element models used in the simulations of delamination 

growth. 


ANALYSIS 


Results from Fatigue Onset and Growth Analysis 

In Figures 8 to 14, the delamination length, a, is plotted versus the number of 
cycles, N, for different input parameters and models. For all results shown, the 
analysis stopped when a 10,000,000 cycles limit - used as input to terminate the 
analysis - was reached. 

Initial Results 

Initial results, as plotted in Figure 8, were obtained using the specified default 
values as input parameters (see appendix of reference 6). For better visualization of the 
results and to be able to identify the differences in the results, the scale on the vertical 
axis was expanded. The results obtained from two-dimensional plane strain (open blue 
circles and solid blue line) and plane stress (open red squares and solid red line) 
models as well as full three-dimensional models (open green diamonds and solid green 
line) were within 4% of the benchmark results (grey crosses and solid grey line). For 
all results shown, the predicted onset occurs prior to the benchmark result onset, 
Nd= 150 cycles. The onset value is lowest for the plane stress results, followed by the 
results obtained from 3D solid models and the plane strain results. This sequence can 
be explained by the computed energy release rate, where the values obtained from the 
plane stress model are slightly higher compared to the results obtained from 3D solid 
and plane strain models as discussed in reference 6. During the growth phase of the 
analysis, the results lie on curves with nearly the same slope parallel to the benchmark 
which suggests that the Paris Law was implemented correctly and is - as expected - 
independent of the model. For all models, the threshold cutoff, where delamination 
growth is terminated and the delamination length remains constant, is predicted close 
to the number of cycles defined by the benchmark. 

Variation of Input Parameters 

Input parameters were varied to study the effect on the computed onset and growth 
behavior during the analysis and results are plotted in Figures 9 and 10. For this 
parametric study, only models with a refined mesh (Aa=0.25 mm) made of plane strain 
elements (CPE4) were used. Once a set of parameters was established that yielded 
good results, the effects of mesh size and element type on the results were studied. 

First, the effect of the initial time increment used in the analysis settings was 
studied as shown in Figure 9 (see appendix of reference 6 for details). The initial time 
increment was varied between /'o=0.01 (one tenth of a single loading cycle, 4=0. Is, 
open blue circles and solid blue line) and /'o=0.000 1 (one thousandth of a single 
loading cycle, open green diamonds and solid green line). For larger initial time 
increments, the onset of delamination shifted towards a lower number of cycles. 
Reducing the initial time increments, however, significantly increased the 
computation time. Based on the results, it was therefore decided to use an initial time 
increment of /'o=0.001 (open red squares and solid red line) for the remainder of the 
study to save computation time. This step is justified by the fact that the results 



obtained for zo=0.001 were almost identical to the values obtained from the analysis 
where a smaller initial time increment was used (z‘o=0.0001). 
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Figure 8. Computed delamination onset and growth: Detail of initial results. 
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Figure 9. Computed de lamination onset and growth obtained for different initial 

time increments. 


Second, the input to define the cyclic loading was varied. In order to define the 
cyclic applied displacement, <5/2, a set of parameters need to be determined as shown 
in the equation (7) and highlighted in Figure 10. Selecting Aj=0 results in a sine curve 
representation of the cyclic load (analysis results are shown in Figure 10). Selecting 
Bf= 0 results in a cosine curve representation of the cyclic load. The selection of the 
starting time, to, causes a phase-shift. Further details about the input parameters are 
discussed in detail in the appendix of reference 6. As shown in Figure 10, the number 
of cycles to delamination onset is the most sensitive to input variations. The results 
obtained for the growth phase however, lie on curves with nearly the same slope, 
parallel to the benchmark. Also the threshold cutoff, where delamination growth is 
expected to stop and the delamination length remains constant, is predicted close to 
the number of cycles defined by the benchmark. A narrow band of results (green solid 
lines), was obtained when the number of terms used to define the Fourier series was 
increased to 50. A Fourier series is used in ABAQUS® Standard during the analysis 
to approximate the periodic cyclic loading (see appendix). The narrow band of 
results was in good agreement with the benchmark curve (grey solid line) compared to 
the results obtained for the default setting of 1 1 Fourier terms which showed more 
scatter (black dashed lines). It should be noted that one term in a Fourier series is 
sufficient to exactly represent a simple periodic function such as the simple sine 
function used in the current example. At the same time, a large number of terms 
should improve the approximation of a more complicated cyclic load as demonstrated. 
Good results could only be obtained when the analyses were performed with 
ABAQUS 8 Standard 6.9EF. Earlier versions showed unexplainable wide variations in 
results. The study was repeated for Bi = 0 which yielded similar results as discussed in 
detail in reference 6. 
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Figure 10. Computed delamination onset and growth obtained for sine representation 

of the cyclic loading. 


Based on the current results, it was decided to use a sine curve representation 
(A i=0) in combination with the starting time, G=0.0 for the remainder of the study. 

Third, the release tolerance was varied. Once a user specified release tolerance 
((G -G c )/G c > release tolerance) is exceeded during the analysis in ABAQUS' 
Standard, a cutback operation is performed which reduces the time increment. The 
cutback reduces the degree of overshoot and improves the accuracy of the local 
solution [10]. In the current study, varying the input between 0.2 (the default value) 
and 0.01 did not have any effect on the computed onset and growth behavior during 
the analysis. It is assumed that the release tolerance only affects static delamination 
propagation and not growth under cyclic loading studied here. 

Fourth, the solution controls were varied. It is generally not required nor 
recommended to modify the solution controls in ABAQUS 8 Standard. Since some of 
the solution controls, however, were modified in the DCB example problem provided 
with the ABAQUS* installation, it was decided to briefly address the effect on the 
computed onset and growth behavior during the analysis. The first input parameter 
that defines the iteration number at which the periodicity condition is first imposed, 
was not modified and the default value was kept. A series of analyses were performed 
where the other four input parameters were varied between 100 and 10" 4 . Changing the 
input did not have any effect on the computed onset and growth behavior during the 
analysis, however it significantly influenced the computation time. The analysis 
required only about a tenth of the computation time when the second input parameter 
was modified from the default value (5T0' 3 ) to values larger than 10" 1 . Changing the 
other remaining three parameters did not have any effect for this example. Further 
details about the input parameters are discussed in reference 6, where also sample 
input files are provided. 

Variation of Mesh Size and Element Type 

The results obtained for models with different mesh sizes and different element 
types are shown in Figures 11 to 13. For models made of plane strain elements 
(CPE4), the element length, Aa, at the delamination tip was varied as discussed earlier. 
The results obtained from the respective models are plotted in Figure 11. Excellent 
agreement with the benchmark curve (grey crosses and grey solid line) could be 
achieved for element lengths up to Aa= 1.25 mm, as shown in Figure 11. For element 
length Aa=\.61 mm (purple triangles and solid purple line), the predicted onset occurs 
for a slightly higher number of cycles. The observed mismatch is largely due to the 
increased element length, which causes the first growth step to be larger. The 
following growth increments are relatively large due to the increased element length 
(Aa=\.61 mm). However, during the growth phase of the analysis, the results from all 
models he on curves with the same slope parallel to the benchmark, which suggests 
that the Paris Law was implemented correctly and is, as expected independent, of the 
mesh. For all models, the threshold cutoff, where delamination growth is expected to 
stop and the delamination length remains constant, is predicted close to the number of 
cycles defined by the benchmark. The total computation time was between 70 s for the 
coarsest mesh and 1030 s for the finest mesh as shown in Figure 12 1 . 


1 CPU time on Dual-Core AMD Opteron(tm) Processor 8220 SE 
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Figure 11. Computed delamination onset and growth behavior for different plane 

strain models. 
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Figure 12. Required analysis time for models with different crack tip element lengths. 


The analyses were repeated for models made of plane stress elements (CPS4), 
where the element length, Aa, at the delamination tip was varied as before. The results 
obtained from the respective models are discussed in reference 6. The computation 
times are included in Figure 12 for comparison. 

The results obtained for models made of 3D solid and continuum shell elements 
showed the same trends and are discussed in detail in reference 6. The computational 
effort, however, was reduced by a factor of 2.5 to 2.9 for analyses performed using the 
continuum shell element models. For the models made of 3D solid elements, the total 
computation time was between 5800 s for the coarsest mesh and about 7 days for the 
finest mesh. For the corresponding models made of continuum shell elements, the total 
computation time was between 2000 s for the coarsest mesh and about 2.8 days for the 
finest mesh 1 . The computation times are also included in Figure 12 for comparison. 

For final comparison, finite element analyses were repeated with two-dimensional 
and three-dimensional models with the same element length, Aa= 0.5 mm, at the 
delamination tip. The results obtained from two-dimensional plane strain (open blue 
circles and solid blue line) and plane stress (open red squares and solid red line) 
models as well as full 3D solid models (open green diamonds and solid green line) and 
continuum shell elements (orange x’s and solid orange line) were within 1% of the 
benchmark curve (grey crosses and solid grey line) as shown in Figure 13. The results 
obtained from the continuum shell element model (orange x’s and solid orange line) 
were almost identical compared to results obtained from the full 3D solid model (open 
green diamonds and solid green line) as mentioned above. The results obtained from 
plane stress models (open red squares and solid red line) are close to the results 
obtained from three-dimensional models. 
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Figure 13. Computed delamination onset and growth behavior for different 2D and 3D 

models. 


For all results shown, the predicted onset occurs first for the plane stress results 
and last for the plane strain results. Also, the delamination length obtained from plane 
strain models are slightly lower. This observation can be explained by looking at the 
computed energy release rate. The plane stress model yields a slightly higher energy 
release rate compared to the 3D solid and the plane strain model. Details are discussed 
in reference 6. Delamination onset would therefore occur first in plane stress models 
as indicated by the results plotted in Figure 13. During the stable growth phase, 
however, the results for all models lie on curves with nearly the same slope parallel to 
the benchmark as mentioned earlier. For all models, the threshold cutoff, where 
delamination growth is assumed to stop and the delamination length remains constant, 
is predicted close to the number of cycles defined by the benchmark. 


SUMMARY AND CONCLUSIONS 

The development of a benchmark example for cyclic delamination growth 
prediction is presented and demonstrated for the commercial finite element code 
ABAQUS® Standard. The example is based on a finite element model of a Double 
Cantilever Beam (DCB) specimen, which is independent of the analysis software used 
and allows the assessment of the delamination growth prediction capabilities in 
commercial finite element codes. First, the development of a benchmark example for 
delamination fatigue growth prediction was presented step by step. The number of 
cycles to delamination onset was calculated from the material data for mode I fatigue 
delamination growth onset. The number of cycles during stable delamination growth 
was obtained incrementally from the material data for mode I fatigue delamination 
propagation. For the combined benchmark case of delamination onset and growth, the 
delamination length was calculated for an increasing total number of load cycles. 
Second, starting from an initially straight front, the delamination was allowed to grow 
under cyclic loading. The number of cycles to delamination onset and the number of 
cycles during stable delamination growth for each growth increment were obtained 
from the analysis. 

The results showed the following: 

• In general, good agreement between the results obtained from the growth analysis 
and the benchmark results could be achieved by selecting the appropriate input 
parameters. However, selecting the appropriate input parameters was not 
straightforward and often required an iterative procedure. 

• The onset prediction appeared much more sensitive to the input parameters than 
the growth prediction. 

• Consistent results were obtained when input parameters were selected such that 
50 terms in the Fourier series were used during the execution of ABAQUS 8 
Standard to approximate the periodic cyclic loading. 

• Good agreement between analysis results and the benchmark could be achieved 
when the initial time increment used in the analysis was about one tenth of a 
single loading cycle. 

• Best results were obtained when a sine curve representation of the cyclic applied 
displacement was selected in combination with the starting time, to=0.0. 

• The release tolerance did not have an effect on the analysis or the results. 



• Accurately computing the onset and growth required fine meshes with an 
element length at the tip Aa< 1.0 mm. 

• The solution controls in ABAQUS R Standard had to be modified in order to 
reduce computation time. Even with carefully selected input parameters, the 
analyses for three-dimensional models of a simple DCB specimen required days. 

• Although implemented in ABAQUS® Standard 6.8, version 6.9EF was required to 
obtain consistently good results. 

• Improvements are needed to make this analysis applicable to real case scenarios 
such as more complex specimens or structural components. 

Overall, the results are promising. In a real case scenario, however, where the 
results are unknown, obtaining the right solution will remain challenging. Further 
studies are required which should include the assessment of the propagation 
capabilities in more complex mixed-mode specimens and on a structural level. 

Assessing the implementation in one particular finite element code illustrated the 
value of establishing benchmark solutions since each code requires specific input 
parameters unique to its implementation. Once the parameters have been identified, 
they may then be used with confidence to model delamination growth for more 
complex configurations. 
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